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Abstract 

A new computer code has been developed for predicting the turbulent, and chemically react- 
ing flows with sprays occurring inside of a stratified-charge rotary engine (SCRE). The solution 
procedure is based on an Eulerian-Lagrangian approach where the unsteady, three-dimensional 
Navier-Stokes equations for a perfect gas-mixture with variable properties are solved in generalized, 
Eulerian coordinates on a moving grid by making use of an implicit finite-volume, Steger- Warming 
flux vector splitting scheme, and the liquid-phase equations are solved in Lagrangian coordinates. 
Both the details of the numerical algorithm and the finite-difference predictions of the combus- 
tor flowfield during the opening of exhaust and/or intake, and also during fuel vaporization and 
combustion, are presented. 


I. Introduction 

The rotary combustion engine (RCE) would be desirable as a powerplant for light aircraft, 
drones (including high-altitude application), auxiliary and ground power units, and also for ma- 
rine and industrial application, if only its efficiency could be improved closer to that of diesel 
engines. It has inherent advantages over reciprocating engines in terms of higher airflow capacity, 
higher power-to-weight ratio, and less vibration, among others. An initial attempt to introduce 
a gasoline-fueled rotary engine into the general aviation market was unsuccessful because of poor 
fuel economy, uncertain availability of avgas, and marginal weight advantage over contemporary 
reciprocating engines. Subsequent research sponsored by industry, NASA, and the Navy has led to 
the development of the stratified-charge rotary engine (SCRE) concept, which has demonstrated 
multifuel capability. Current R&D sponsored by NASA is aimed at reducing the cruise brake spe- 
cific fuel consumption (BSFC, lb/bhp-hr) from a current value of about 0.42 to 0.35 or less by the 
end of 1992. The expected improvement will be enabled by further CFD - driven fuel injection, 
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spray, and nozzle optimizations, rotor pocket and nozzle relocations, and related modifications. In 
this paper, we address the CFD aspect of the technology development program for predicting the 
complex flow patterns occurring inside of a Wankel engine. 

Early modelling efforts on the Wankel engine were based on thermodynamic models 1,2 and also 
on one-dimensional modelling of premixed-charge combustion. 3 Multi-dimensional models of the 
Wankel engine are very recent in origin; Grasso et al. 4 have presented the first three-dimensional 
computations of a SCRE during the early stages of flame propagation. Subsequent computations 
performed by Abraham and Bracco 5,6 have led to some important design changes in the rotary 
engine development at John Deere k Co., especially in the fuel injector configuration. Their 
code, REC-3D-FSC-86, is a modifed version of the KIVA code developed at Los Alamos National 
Laboratory 7 for the modelling of reciprocating engines. KIVA makes use of a conditionally stable 
algorithm, and the stability of the KIVA scheme is improved by making use of an acoustic subcycling 
step in order to alleviate the stiffness problems arising from compressibility effects. There appears to 
be considerable room for improvement in the code, since it neglects the spatial gradients whenever 
the grid spacing becomes smaller than some predefined value and also requires excessive CPU time 
when the engine speed becomes small. Shih et al. 8 presented the first two-dimensional computations 
of a motored Wankel engine in the absence of combustion. Their code, LEWIS-2D, is based on 
the Beam- Warming type of ADI method. Their computations have subsequently been extended to 
three dimensions in Steinthorsson et al. 9 Linear stability analysis has shown that the ADI method 
is unconditionally stable in two dimensions but is unconditionally unstable in three dimensions. 
Although artificial dissipation has some stabilizing effect, an excessive amount can impair stability 
and reduce accuracy and convergence. Recently, Li et al. 10 have modified their LEWIS-3D code 
based on upwind schemes together with the incorporation of a k — € turbulence model. 

The present solution procedure differs from both REC-3D-FSC-86 and LEWIS-3D in many ways 
in terms of the numerics, and also the submodels used for turbulence, combustion, and sprays. 


IL Physical Description 

A schematic of the Wankel engine is shown in Fig. 1. The Wankel engine is composed of a 
peripheral housing with provisions for the intake and exhaust ports, fuel injector and spark igniter, 
a three-flank rotor, and a crank shaft. The contour of the inner surface of the outer casing of the 
Wankel engine is composed of a two-lobe peritrochoid. 11 The contour of a rotor revolving along 
an outer housing is represented by a peritrochoid inner envelope. The geometric analysis of the 
rotor and housing surfaces can be found in Yamamoto. 11 The rotor surface is further modified by 
the formation of a rotor pocket. The presence of a rotor pocket not only alters the expansion 
and/or compression ratio of the engine but also plays an important role in modifying the flowfield, 
and mixing and combustion characteristics of the combustor. The present rotor configuration is 
adopted from Steinthorsson et al. 9 
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The rotor turns eccentrically at one third of the crank shaft speed. The three combustion 
chambers of the Wankel engine are the three regions enclosed between the three rotor faces and 
the peritrochoid housing, two side housings, two side seals, and lead and trail apex seals. In 
the present calculations, only one of the three combustion chambers is considered, since leakage 
through the seals is assumed to be negligible. As the rotor revolves around the crank shaft, each of 
the combustion chambers is continually deformed. This produces the necessary compression and 
expansion of the fluid for the required engine performance during each one of the cyclic operations. 
The intake and exhaust ports, spark igniter, and the fuel injector are located along the peritrochoid 
housing as shown in Fig. 1. 

In the present study, the computations are initiated before the opening of the exhaust port 
for the combustor formed with the second rotor flank as shown in Fig. 1. The initial conditions 
correspond to the conditions of quiescent air at pressure, P{ n = 1 atm, and temperature, T* n 
= 300/f. As the rotor moves in the clockwise direction the exhaust port opens and the residual 
gas moves out of the combustion chamber, since the normally imposed pressure in the exhaust 
remains lower than the interior engine pressure during most of the compression cycle. The exhaust 
conditions are given by 


dp __ dyj _ de = 
dn dn dn * 


v = -C dc 


P = Ptzh, U = W = 0, 
2(P-P„»)1 0li 


I 


( 1 ) 


where p> e, and y are the fluid density, internal energy, and mass fraction, respectively; u, v, and w 
are the velocity components in Cartesian coordinates; (= 0.9) is the discharge coefficient; and 
subscript i and n represent species and the normal component of the boundaries, respectively. 

As the rotor moves further in the clockwise direction, the intake port opens and fresh air moves 
into the combustion chamber. There is an overlapping region during which both the exhaust and 
intake ports are simultaneously open before the exhaust port closes completely. During this process 
not only the residual gas but also some of the fresh intake might escape through the exhaust port. 
Most of the intake occurs during the expansion stroke of the engine. The inport conditions are 
given by 


P — Pint y 'P “ Pint* Vi — 

P = Piniy U = W = 0, 
[2 (i^nt — P)1°' 8 


V = -Cdc 
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As the rotor turns further, the intake port closes and the liquid fuel is injected into the chamber 
during the compression stroke before the top-dead center (TDC) is reached. Spark injection pro- 
vides the initial energy needed for the early droplet evaporation and also helps to promote ignition 
of the vaporized fuel and air mixture. Most of the combustion is completed as the rotor moves past 
the TDC, and most of the residual combustion products are eventually driven out of the combustion 
chamber through the exhaust port. The whole combustion performance is determined by a very 
complex interaction of various engine parameters including the location of the exhaust and intake 
ports, shape of the rotor pocket, injector and spark timings, fuel properties, and many others. 

III. Gas-Phase Equations in Generalized Coordinates 


The governing unsteady equations based on the conservation of mass, momentum, energy, and 
species for turbulent, reacting, and compressible flows are presented in strong conservation law form. 
The exchanges of mass, momentum, and energy through liquid-phase interaction are considered by 
the inclusion of appropriate source terms. The Reynolds- averaged equations are formulated in 
generalized coordinates to accommodate the time-variation of the complex combustor geometry. 
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and x , y, and z are the Cartesian coordinates in the physical space; r) } and f are the coordinates 
in the computational space; D is the determinant of the matrix, J in Eq. ( 6 ), and is also a 
measure of the volume of a computational cell; y, is the mass fraction of the ith species; S c is 
a vector representing the source terms arising from the finite-rate chemical reactions; W+ is the 
molecular weight of the species; 14 is the stoichiometric ratio of the ith species participating in 
a given reaction step; A and E a are the pre-exponential coefficient and activation energy of a 
given Arrhenius reaction-rate term; 5/ is a vector representing the source terms arising from the 
liquid-phase interaction; n* is the number of droplets in a fcth characteristic representing a group 
of droplets; m * is the vaporization rate of a droplet belonging to the kth characteristic; r* is the 
droplet radius; hj 9 and lk,eff are the enthalpy of the fuel vapour at the droplet surface, and the 
effective latent heat of vaporization; fit is the turbulent viscosity; ki m and fn m are the thermal 
conductivity and laminar viscosity of the gas mixture and are determined using Wilke’s mixing 
rule with fourth-order polynomial fits based upon temperature dependence 12 ; C pm is the specific 
heat of the gas mixture at constant pressure and is also determined from fourth-order polynomial 
fits involving temperature dependence; Pr* (= 0.90) is the turbulent Prandtl number; Set (= 0.90) 
is the turbulent Schmidt number; the subscripts /, o, /, l y c, m, and k represent fuel, oxidizer, 
liquid-phase, laminar, chemical reaction, gaseous mixture, and characteristic, respectively. 

The pressure and temperature are calculated iteratively from the following procedure: 
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where h j- is the heat of formation of ith species, and J2 U is the universal gas constant. Equation 

(5) is the equation of state for a gas mixture of N t species. The Jacobians of the coordinate 

transformation are given by 
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The metric coefficients resulting from the coordinate transformation are evaluated from the 
following identities: 
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It is noteworthy that tfye following equations represent the metric invariant terms arising from the 
coordinate transformation: 
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When the governing equations are formulated in strong conservation form, it is essential that 
the left-hand side of Eqs. (9) to (12) vanish identically when the derivatives are approximated by 
finite-differences; otherwise spurious source terms may result from geometrically induced errors. 13 
Equations (10) to (12) are satisfied identically when central differences are used to evaluate the 
spatial derivatives. This is true since the metric identities in Eq. (8) are written in conservative 
form. However, the determinant of the coordinate transformation is computed numerically from 
the solution of Eq. (9) in order to avoid grid-motion induced errors. 13 

IV. Details of the Spray, Combustion, and Turbulence Models 


Here we provide a brief description of the spray model as it is adopted from Raju and 
Sirignano. 14,15 The solution of the liquid-phase equations is extended further in the present study 
from the two-dimensional to three-dimensional computations. Also, several modifications are incor- 
porated into the interpolation procedure between the Eulerian and Lagrangian coordinate systems 
as the gas-phase computations are performed in the generalized coordinates as opposed to the 
Cartesian coordinates used in Raju et al. 14,16 The interaction between the two phases is taken 
into account by the following procedure: (1) In order to obtain the solution of the liquid-phase 
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equations, it is first necessary to know the gas-phase properties at the particle locations. In the 
present computations, the gas-phase properties are evaluated by using a second-order accurate in- 
terpolation method involving volume-weighted averaging; (2) The ordinary differential equations 
describing particle size, position, and velocity are solved by the second-order accurate Runge-Kutta 
method. The partial differential equation describing the transient temperature variation within the 
droplet interior is based on a simplified vortex model and is solved by an implicit method. The 
formulation for the droplet vaporization rate is based either on a simplified gas-phase boundary 
layer analysis or on a simplified correlation, 15 depending upon the droplet Reynolds number; (3) 
Finally, after the liquid-phase equations are solved, the source terms evaluated at the particle lo- 
cation are redistributed amongst the eight computational nodes surrounding the particle by using 
volume-weighted averaging. 

The success of the spray model also depends a great deal on the correct specification of the 
injector exit conditions. The injector is an eight-holes configuration and is located along the middle 
of the peritrochoid housing as shown in Fig. 1. The timing of the injector opening and closing is 
determined by the given engine operating conditions. The fuel emerges in a fan shape consisting of 
eight streams as shown in Figs. 5 and 6. Both the initial droplet velocities and temperatures are 
assumed to be known, and the droplet sizes are determined by the Rosen-Ramler distribution. 16 
The droplet injection timing is determined by the resolution of the computational cells used in 
the gas-phase computations. 14 The present model does not take into account details of the liquid 
filament breakup and its subsequent effect on the conditions at the injector exit. 

The solution procedure could perhaps be improved with the consideration of droplet dispersion 
due to turbulence. However, the effects of turbulent dispersion in the modelling of combusting 
sprays were found to be small in a previous study, in comparison to the uncertainties in the 
specification of the initial conditions at the injector exit. 17 The present model is based on a dilute 
spray approximation where the spray characteristics are based on an isolated droplet behaviour. 
O’Rourke and Bracco, 18 Greenberg and Tambour, 19 and Asheim et al. 20 have modelled liquid sprays 
including droplet collisions. The importance of droplet collision and breakup in the overall spray 
behaviour is not well established, especially in the regions of the spray where the droplet loading 
is low. In the present computations, the effect of variable properties in the liquid-phase is not 
considered, though this factor becomes very important when the droplets vaporize near the critical 
conditions. 21,22 

The combustion model is based on an analogous treatment of laminar diffusion flames with the 
assumption that no envelope flame exists. The reaction rate is determined based on a single-global 
kinetic mechanism of Westbrook and Dryer. 23 For n-decane, the kinetic mechanism is given by 

C 10 H 22 + 15.5(02 + 3.76# 2 ) -+ 10 C0 2 + 11 H 2 0 + 58.28^2 (13) 

By assuming equal binary diffusivities for all the species in the mixture, the concentrations of 
N 2 > C0 2 , and H 2 0 can be determined from simple algebraic relationships based on the atomic 
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balance of the constituent species, once the mass fractions of fuel and oxidizer are known from the 
solution of the two gas-phase equations based on the conservation of fuel and oxidizer. 

Vh 2 o = Ki - K\Kiyo 2 - K2yc l0 H 22 

yco 2 = K 2 KZ - K\KiK$yo 2 - K^Kzyc X0 H 22 ( 14 ) 

Vn 2 — 1 ~ - yo 2 (l - K\Ki - K\KiK$ - ycw^aU ~ K* " -KjAs) 

where K\ = 4.29, If* = 0.08723, and K$ = 2.222. Note that the effect of turbulence on the 
reaction rate can be very important, but is not considered in the present solution procedure since 
a realistic model is not currently available. In the near future, we are planning to implement the 
turbulence-reaction model used in Raju and Sirignano 14,15 which is based on the eddy break-up 
model of Spalding 24 requiring the solution of an additional equation involving the square of fuel 
concentration fluctuations. The eddy break-up model provided some useful results in the modelling 
of premixed flames, however, its applicability in a spray environment is uncertain. 

The turbulence model used is a constant eddy viscosity model of Steinthorson et al. 9 where the 
turbulent diffusivity is given by 


\i t = a T Clp (15) 

and ax is a function of the crank angle 0, and fi is the crank speed. One obvious discrepancy of 
this model is its failure to satisfy the condition of /it — 0 at the walls. 9 The k — e turbulence model 
of Launder and Spalding 25 will soon be incorporated into our solution procedure. 

V. Details of Flux Vector Splitting 


The present finite-difference formulation is based on an upwind scheme because of its superior 
numerical stability, and efficiency properties compared to those of a centered difference scheme. 26 
The most widely used flux vector splitting methods are those of Steger and Warming, 26 van Leer, 27 
and Roe. 28 Recently, there is a considerable interest in extending these methods for the modelling of 
reactive flows to solve problems emerging from the design of the National Aerospace Plane (NASP), 
and Air-assisted Orbital Transfer vehicles (AOTV). 29,30 Details of the generalization of the Steger 
and Warming flux vector splitting for a perfect gas mixture with variable properties are presented 
in this section. Since the flux vector F(<2) of Eq. (1) retains its homogenous property for the 
equation of state considered, the flux vector can be split into two parts, 

F = + F~ (16) 

where is the subvector associated with the non-negative eigenvalues of A , F“ is the subvector 
associated with the non-positive eigenvalues of A, and A is the Jacobian matrix, |^.The eigenvalues 
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of the matrix, A are given by 
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and 


h = ho+ r C po dT - £ W 0 C pm T - ^(u 2 + v 2 + u; 2 ), 
J T ref Hu * 

h = hi+ ^ C p idT - ^-W x C pm T, 
hi = h 2 + f T C p2 dT - -£-W 0 C pm T, 

JT ref K u 

and also *y — Cp m j C vmy = = and U ”b 

The resulting components of the split fluxes F ± are given by 
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F f = [*f ^ (u + La) + ^(« “ £«<*)] 

Ff = %^P D [ A l 3^ w + ^( v + &•) + ( v ““ M 

- ± - ± - ( 20 ) 

^ 4± = [ A f^r w + ^( w + £*°) + nM w - f«°)] 

if = fo^[Af^(fc-C^.T) + ^(fc + ^a) + ^(fc-&a)] 

if = Vi*? 

if = »*f 

In this section, we have presented the derivation of the the split fluxes associated with the flux 
vector F(Q). The corresponding split fluxes associated with the vectors G(Q) and H(Q) can be 
derived in a similar way. 

VI. Details of the Numerical Method 

Solution for the gas-phase equations is obtained by making use of a finite-volume, Lower- 
Upper (LU) decomposition scheme. The governing equations are linearized in a delta form where 
the nonlinear terms associated with finite-rate chemistry and convection are treated implicitly, 
while the diffusion terms and the source terms arising through liquid-phase interaction are treated 
explicitly. The time-linearized governing equations in delta form are 

[D n+1 I + At (Sf A + + S~ B + + 6~C + + 8} A~ + S+ B~ + S+C~ - L)]AQ 

= -Q n AD + Atr ( 21 ) 

where A ± = Sjjr, B ± = C ± = L= ^4, At is the time step size, 6 + and 6~ are forward 
and backward differences, respectively, and 

AQ = Q n AD + D n+1 AQ, 

r = -Sf F + - S~G + - 5~ H + - 6{F~ - 8+G~ 

-6+H- + 6tF v + Sr,G v + 6 ( H v + S e + S t ( 22 ) 
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(23) 

(24) 


Upon factoring Eq. (20) we obtain the following sequence: 

[£>" +1 / + At (5 £ + A _ + 8+B~ + 6+C~ - L)] A Q* = -Q n AD + Atr 
[i) n+1 / + At (5f A + + 6~B + + S~C + )] A Q = I D n+1 A Q* 

It is noteworthy that to be consistent with the objective of deriving a finite-volume code, 
the split-flux differences in Eq. (21) are implemented according to Monotone Upstream-Centered 
Schemes for Conservation Laws (MUSCL)-type differencing. 31,32 The fluxes at the cell faces are 
first obtained by a fully upwind first-order accurate interpolation, and then centered differences are 
used for both the forward and backward spatial operators evaluated at the cell centers. Centered 
differences are also used for evaluating the spatial operators associated with the viscous terms. 

For the dynamic grid calculations, the metric quantities are evaluated at time level n+1, and 
£jn+i evaluated from the solution of Eq. (9) by using an explicit method. The numerical grid 
is generated by an algebraic technique 8 with the help of the grid-generation code taken from the 
LEWIS-3D code. 9 

By adopting an algorithm taken from the RPLUS-3D code, 12 the present code is vectorized 
rather efficiently by operating on all points in a diagonal plane of the computational space, simul- 
taneously. The diagonal plane is one on which i+j+k = constant. The integration proceeds during 
the backward and forward substitution steps from one corner node of the computational grid and 
ends at a corner node which is farthest from the initial corner node. 

The boundary conditions are implemented explicitly by defining a layer of phantom cells outside 
the boundaries of the computational domain. It is also noteworthy that the left-hand side operators 
of Eqs. (23) to (24) require block diagonal inversions. 

VII. Results and Discussion 

Here we present the results of our preliminary computations for a single case corresponding 
to the operating conditions listed in Table 1, where subscripts o, c, r, and h represent opening, 
closing, rotor, and housing, respectively, is the initial droplet velocity, and T^i is the initial 
droplet temperature. Isothermal wall conditions are implemented in the case considered. 
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Table 1. Operating Conditions 

Engine Parameters 
(see Fig. 1) 

Generating Radius(R) = 0.1064 m 
Eccentricity (E) = 0.01542 m 
Clearance(C) = 0.004 m 
Chamber Width(W)= 0.07 m 
Port Width(W P )= 0.05 m 

Engine Speed 

4000 rpm 

Intake Port 

0 o = -1.26 rad, 9c = 5.96 rad, Pint = 1.25 atm, 
Tint = 300 K, Y f ,int = 0 

Exhaust Port 

$o = -5.96 rad, 0 C = 1.07 rad, P ex h = 0.85 atm 

Fuel Injector 

$o = 8.3 rad, 9c — 8.75 rad, V^,- = 100 m/s, 
Tn = 300 K 

Spark ignition Timings 

0 o = 8.35 rad, 6c = 8.475 rad 

Temperature of Rotor 
and Housing Surfaces 

T h = T r = 400 K 


The computations are performed with a variable time-step corresponding to a maximum CFL 
number of 20, and on a grid with a mesh size of i=31, j= 16, and fc=20, where *, j, and k represent 
the coordinate surfaces in the direction extending from the trailing-edge surface to the leading-edge 
surface of the combustor, from the rotor to housing surface, and from the side wall to the symmetry 
plane of the domain between the end-to-end side walls, respectively. 

As described in Section II, the computations are initiated before the opening of the exhaust port, 
and are terminated after the completion of combustion process. Figure 2 shows the variation of 
the engine volume versus crank angle. The indicated engine volume is obtained from the numerical 
integration of the individual computational-cell volumes. The figure demonstrates the ability of 
the numerical method to accurately reproduce the combustor volume changes corresponding to 
a maximum and minimum displacement volume of 750 c.c. and 115 c.c, respectively, yielding a 
compression ratio of about 6.5. 

Figure 3 shows flow patterns (particle traces) during the middle of the exhaust (Fig. 3a), of 
the intake (Fig. 3b), and of the simultaneous opening of both exhaust and intake ports (Fig. 3c) 
together with the schematic of a Wankel engine (Fig. 3d). The particle traces are coloured based 
upon the local value of the internal energy e. Fig. 3a shows clearly that the bulk fluid motion of the 
residual gas is essentially directed out of the engine chamber through the exhaust port. The fluid 
motion is clearly seen to be influenced by the combined effect of the positive pressure difference 
that exists between the chamber and back pressure, and of the rotor motion on the fluid near the 
rotor surface. Figure 3c shows that while part of the residual gas is escaping through the exhaust 
port, fresh air is entering the chamber through the intake port. Because of the difference that exists 
between the intake and exhaust pressure (1,25 and 0.85 atm), part of the fresh air emerging from 
the intake is drawn towards the exhaust before it recirculates in the rotor pocket. 
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Figure 3b shows a complex flow pattern created by a strong jet of fresh air emerging from the 
intake in a cross flow. The cross flow is created by the clockwise movement of the rotor surface. 
The flow pattern reveals the existence of two clearly defined recirculation (low-pressure) regions on 
both sides of the intake along the tth direction. Upon the impingement of the intake jet on the 
rotor surface, the jet, near the leading edge, is partly drawn into the low-pressure region, and part 
of the fluid is carried over backwards through the openings between the outer edge of the jet and 
the side- walls. 

Figure 4, which is similar to Fig. 3b, shows an experimental flow-visualization result obtained 
by Hamady et al. 33 using a transparent-sided rotary engine motoring test rig. “Microballoon” 
seeding material was illuminated by a pulsed laser light and recorded by high-speed photography. 
At the leading edge of the jet, a recirculatory flow pattern quite similar to that shown in Fig. 
3b is clearly seen. Although quantitative comparison has not yet been attempted, the degree of 
qualitative agreement noted here is quite encouraging. The trailing vortex, which is clearly shown 
in Fig. 3b, is beyond the field of view of Fig. 4. 

The droplet trajectories at 6 - 8.5 and 8.75 rad are shown in Figs. 5 and 6. The polydisperse 
character of the spray is represented by different sized circles which are indicative of the size of 
the initial droplets. The initial droplet sizes range between 10 pm < < 30 /im, and the initial 

droplet Reynolds numbers vary between 75 < Re^i < 600. The wide disparity in Re^i is a result of 
the steep rise in the chamber pressure due to combustion from 6 to 35 atm during the time of fuel 
injection. Because of the large initial momentum associated with these particles, they retain their 
initial path as described by their initial conditions. It takes about 1.5 ms for the largest particles 
to vaporize. The deflection of the particles in the direction of the gaseous flow is evident from Fig. 
6 as the droplets become smaller due to evaporation. 

The temperature distribution within the combustion chamber during the early stages of flame 
propagation at 6 = 8.5 rad, and after combustion at 6 = 10 rad is shown in Figs. 7 and 8, 
respectively. In Fig. 7, the highest temperature region (2900 K) is confined to the region near 
the rotor pocket, where the liquid fuel is injected. The temperatures are lower near the walls and 
in the clearance regions near the leading and trailing apex seals, where the heat transfer to the 
walls is greatest because of the high surface-to-volume ratio. In Fig. 8, the highest temperature 
region extends all the way from the region near the rotor pocket to the region near the leading 
apex seal, while the temperatures are lower in the region near the trailing apex seal. During the 
expansion stroke, the region near the leading apex seal becomes wider, causing a decrease in the 
heat transfer rate to the walls in that region, while heat transfer to that region within the chamber 
interior increases due to convection of the fluid as influenced by the rotor rotation. An opposite 
trend is observed in the region near the trailing seal. 

The gaseous-fuel mass fraction contours at 9 = 8.5 and 8.75 rad are presented in Figs. 9 
and 10. It is noteworthy that in presenting some of the results involving the iso-contour lines of 
fuel concentration in Figs. 9 and 10, and also of pressure in Fig. 13, some of the contour lines 
representing the near-maximum values are shown in dotted lines. Thus, Figs. 9-10 show that the 
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region near the fuel injector location is fuel rich. Diffusion of the fuel concentration with time and 
the influence of convection on the distribution of the fuel concentration more towards the leading 
region is also evident from the comparison of these figures. While the stratified charge gives rise to 
a diffusion flame, a careful examination also reveals that part of the evaporated-fuel and oxidizer 
mixure burns like a premixed flame. This is demonstrated by an absence of fuel concentration 
in the region near the rotor pocket surface where the fuel concentration results otherwise from 
the presence of liquid fuel in that region, as shown in Figs. 6 and 7. It is more likely that the 
combustion characteristics in that region might be influenced by an isolated-combusting droplet 
behaviour. 

Fig. 11 shows the angular variation of the amount of the evaporated-fuel and also the amount of 
reacted-fuel. The results are obtained by integrating the source term contributions of the gas- phase 
equations arising from the production of fuel due to evaporation and also from the consumption 
of reactants due to combustion. It is noteworthy that the total amount of liquid fuel injected is 
determined based on an equivalence ratio of 0.7. The results show that the total time for complete 
vaporization and also combustion is less than 2 ms. The slope of these curves indicates that most of 
the fuel, after it evaporates, reacts quickly with oxidizer to form products. This in turn implies that 
most of the fuel burns in a premixed-flame environment. Since the vaporization rate and ignition- 
delay characteristics of this model were not known in advance, however, the fuel injection and spark 
timings were arbitrarily chosen to be 66° and 63°, respectively, before TDC. These conditions, 
which in retrospect are clearly non-optimal, correspond to very advanced fuel injection and spark 
timings. Under more optimal engine operating conditions, the fuel injection evidently shov.id not 
begin before 30° to 45° from the TDC. Both the vaporization and combustion characteristics might 
be quite different if these timings are chosen according to optimized operating conditions. Future 
work will address the optimization of these timings in terms of the overall combustion behaviour. 
Note also that because of these advanced timings, additional energy has to be supplied for the work 
to be performed during the remainder of the compression process from the time after combustion 
to the time before the TDC is reached. 

Figure 12 shows the velocity vector plots at four different crank angles: Fig. 12a at the beginning 
of fuel injection, Figs. 12b and 12c during combustion, and Fig. 12d after combustion. In these 
velocity vector plots, only three different sizes of arrow symbols are used to distinguish the variation 
between the maximum and minimum values in magnitude. These plots indicate that the direction 
of fluid motion near the symmetry plane is mainly determined by the rotor motion; however, for 
a brief period during the early stages of flame propagation as shown in Fig. 12b, the expanding 
gases do create a motion in a direction opposite to the main bulk flow. The non-uniformity in 
pressure between the leading and trailing regions as shown by the pressure contours in Fig. 13 
better explains the reason for the strong bulk fluid motion created by the rotor movement. Recall 
that, as in Figs. 9 and 10, some of the near-maximum contours are shown in dotted lines. These 
figures do show a more or less gradual decrease of pressure from the trailing to the leading apex 
seals, and the pressure decreases in the direction of the fluid motion. The higher pressure drop, as 
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expected from the substantial friction losses in the regions closer to the apex seals, is also evident 
from these figures. 


VIII. Concluding Remarks 


We have presented a description of a new computer code developed for the modelling of 
stratified-charge rotary engine performance based on the solution of the unsteady, three-dimensional 
Navier-Stokes equations, with the use of convenient submodels for turbulence, combustion, and 
sprays. The details of the rotary engine flowfield during exhaust and/or intake processes and 
compression stroke, and also the details of the mixing, vaporization processes during and after 
combustion have been presented for a single case with advanced fuel injection and spark ignition 
timings. The salient features of this work are summarized below: 

1. The code takes approximately 3 CPU-hours, when the calculations are performed on a grid with 
a mesh size of 31x16x20 on a CRAY Y-MP, for a non-reacting case, and it takes about 7.5 to 10 
CPU-hours for a reactive case with sprays. For the non-reactive case, the solution can be marched 
in time non-iteratively, but, for a reactive case, the solution is obtained by an iterative procedure. 

2. One apparent finding of our study is that vaporization appears to be more rate-controlling than 
mixing during the combustion process, at least in the case that we have studied with advanced fuel 
injection and spark ignition. 

3. There is a good degree of qualitative agreement between the prediction and an experimental 
flow-visualization pattern of Hamady et al. 33 obtained during the intake process. 

4. The present solution procedure makes use of an extremely simplified constant diffusivity turbu- 
lence model. The k — e turbulence model of Launder and Spalding 25 will soon be incorporated into 
our solution procedure. 

5. The present combustion model is based on laminar kinetics. While recognizing the fact that 
no realistic model would be available in the foreseeable future for a proper treatment of the ef- 
fect of turbulence on combustion reaction rates, the eddy break-up model of Spalding 24 will be 
incorporated into our solution procedure in order to account for some effect of turbulence on the 
combustion processes. 

6. After implementing the above-mentioned modifications we will conduct a parametric study in 
order to optimize the location and also timing of the fuel injector and spark igniter. 
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Fig. 1 Schematic of the Wankel engine that was studied. 











ORIGINAL PAGE 

COLOR PHOTOGRAPH 



25 


Flow patterns during exhaust and/or intake processes. 
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Fig. 4 Experimentally obtained flow pattern during intake 
at 1170 rpm in a supercharged rotary engine. 
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DROPLET TRAJECTORIES 
CRANK ANGLE= 8.5 RAD. 


Fig. 5 Droplet trajectories at 6 = 8.5 rad. 
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Temperature distribution at Q — 8.5 rad. 
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Temperature distribution at 6 = 10 rad. 



FUEL CONCENTRATION CONTOURS. 
CRANK ANGLE' 8.5 RAD. K=19 
CMAX' 0 . 1 4 36 . CM I N = 0.0000.DC= .010 





FUEL CONCENTRATION CONTOURS. 
CRANK ANGLE' 8.5 RAO. J» 7 
CMAX= 0.0290.CMIN' 0.0000.DC= .002 


Fig. 9 Fuel mass fraction contours at 6 = 8.5 rad. 
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FUEL CONCENTRATION CONTOURS. 
CRANK ANGLE = 8.7 RAD. K = 19 

CMAXr 0 . 1 400 . CM I N= 0 . 0000 . DC= .010 





FUEL CONCENTRATION CONTOURS. 
CRANK ANGLE= 8.7 RAO. J= 7 
CMAX = 0 . 0355 . CMIN= 0.0000.DC= .002 


Fig. 10 Fuel mass fraction contours at 9 = 8.75 rad 
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VELOCITY VECTOR PLOT . 
CRANK ANGLE= 8.3 RAD. K = 19 
VMAX= 30.40 M/S 






Fig. 12 Velocity vector plots. 
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( b ) PRESSURE CONTOURS. 
CRANK ANGLE= 8.7 RAD K : 
. 36E + 07 , PMIN= .30E + 07.DP: 


19 

. 45E+05 



PMAX 


(c) pressure contours, 

CRANK ANGLE= 10.0 raq « = 
13E+07,PMIN= .11E+07DP 


Fig. 13 Pressure contours. 
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